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ABSTRACT 

We study wakes and gap opening by low mass planets in gaseous protoplanetary disks threaded by 
net vertical magnetic fields which drive magnetohydrodynamical (MHD) turbulence through the mag- 
netorotational instabilty (MRI) , using three dimensional simulations in the unstratified local shearing 
box approximation. The wakes, which are excited by the planets, are damped by shocks similar to 
the wake damping in inviscid hydrodynamic (HD) disks. Angular momentum deposition by shock 
damping opens gaps in both MHD turbulent disks and inviscid HD disks even for low mass planets, 
in contradiction to the "thermal criterion" for gap opening. To test the "viscous criterion" , we com- 
pared gap properties in MRI-turbulent disks to those in viscous HD disks having the same stress, and 
found that the same mass planet opens a significantly deeper and wider gap in net vertical flux MHD 
disks than in viscous HD disks. This difference arises due to the efficient magnetic field transport 
into the gap region in MRI disks, leading to a larger effective a within the gap. Thus, across the 
gap, the Maxwell stress profile is smoother than the gap density profile, and a deeper gap is needed 
for the Maxwell stress gradient to balance the planetary torque density. We also confirmed the large 
excess torque close to the planet in MHD disks, and found that long-lived density features (termed 
zonal flows) produced by the MRI can affect planet migration. The comparison with previous results 
from net toroidal flux/zero flux MHD simulations indicates that the magnetic field geometry plays an 
important role in the gap opening process. Overall, our results suggest that gaps can be commonly 
produced by low mass planets in realistic protoplanetary disks, and caution the use of a constant 
a-viscosity to model gaps in protoplanetary disks. 

Subject headings: accretion disks, stars: formation, stars: pre-main sequence 



1. INTRODUCTION 

I> . During the past decade, exoplanet surveys (e.g. 
On ■ Kepler) have greatly advanced our knowledge of the de- 
' mographics of exoplanets (Howard et al. 2012, Dong & 
Zhu 2013), which in turn sheds light on planet forma- 
tion. In the near future, ALMA will directly probe young 
planets still embedded in gaseous protoplanetary disks. 
, Interpreting such observations will require a better un- 
ff^ ' derstanding of how planets interact with gaseous disks. 
y—^ ' In particular, it is well known that planets can excite 
IV , density waves in disks, which can lead to angular momen- 
. ^ ' tum exchange between planets and disks and cause plan- 
ks ' ets to migrate (Goldreich & Tremaine 1979, Ward 1986, 
'jjj ' Tanaka et al. 2002, Kley & Nelson 2012, and references 
. therein). When these excited waves deposit angular mo- 
mentum back to the disk, gaps will be opened by planets 
in disks (Lin & Papaloizou 1979, 1993, Artymowicz & 
Lubow 1994, Takeuchi et al. 1996, Bryden et al. 1999, 
Kley 1999). Gap opening is important for understanding 
both planet migration and growth. After a gap is opened, 
planet migration slows to the disk viscous timescale (Lin 
& Papaloizou 1986, Nelson et al. 2000), which is normally 
referred to as Type II migration. Slowing down of migra- 
tion is important for the survival of planets on the 10^ 
yrs disk life time. Moreover, after a gap opens accretion 
onto the planet is limited to flow from the circumplan- 
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etary disk (Lubow et al. 1999 and Lubow & D'Angelo 
2006). The circumplanetary disk controls both the final 
mass of the planet and the formation of planetary satel- 
lites (Ward & Canup 2010). A better understanding of 
each of these topics motivates the study of gap opening 
in more realistic models of protoplanetary disks. 

Conventionally, it is thought that two conditions must 
be fulfilled simultaneously for a planet to open a gap 
(Lin k Papaloizou 1993, Bryden et al. 1999). The first 
is the "thermal criterion", which states that the planet's 
Hill radius needs to be larger than the disk scale height 
so that the density wave shocks just as it is excited. In 
other words, the planet mass needs to be higher than the 
so called - local disk thermal mass 



Mth = 



Mq 



1/2 



V5AU 
(1) 

where Mj is the mass of Jupiter, is the disk local sound 
speed, and Rp and Qp are the planet's semi-major axis 
and orbital frequency. Recently, this "thermal criterion" 
for gap opening has been questioned, since even small 
amplitude waves can still steepen and shock in disks af- 
ter traveling some distance (Goodman & Rafikov 2001, 
Rafikov 2002a, Muto et al. 2010, Dong et al. 2011b). This 
makes possible gap opening by planets with masses sig- 
nificantly lower than the thermal mass Mth, see Rafikov 
(2002b), Li et al. (2009), DuffeU & MacFadyen (2012). 

The second condition is the "viscous criterion" , which 
states that the disk viscosity must be low enough for the 
torque from the planet to overcome the viscous torque. 
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Combined with the "thermal criterion" , this requires 
(Lin & Papaloizou 1993) 



40z/ 



(2) 



whore is the kinematic viscosity of the viscous HD 
disks. This "viscous criterion" has been analyzed in de- 
tail by Crida ct al. (2006) through studying the stream 
lines and a refined formulae was suggested. 

An obvious difficulty with the "viscous criterion" is 
that protoplanetary disks are in fact inviscid; angular 
momentum transport and accretion arc thought to be 
associated with MHD turbulence driven by the magnc- 
torotational instability (MRI, Balbus & Hawley 1991, 
1998 for a review). In ideal MHD the MRI has been 
extensively studied using both local (Hawley et al. 1995, 
Stone et al. 1996, Miller & Stone 2000, Davis et al. 2010, 
Guan & Gammic 2011) and global numerical simulations 
(Armitagc 1998, Hawley 2000, 2001, Fromang & Nelson 
2006, Flock et al. 2011, Beckwith et al. 2011). Some as- 
pects of MRI-driven turbulence can be represented by an 
"a" disk model (Shakura & Sunyaev 1973) in a statistical 
sense (Balbus & Hawley 1991,1998, Balbus & Papaloizou 
1999): 



a 
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Pocl 



OiRe + Q-Max 



{{pVxV'y}) {{B,By)) 



(3) 



where Tmri is the averaged xy component of the to- 
tal stress tensor, an,-, is the normalized Reynolds stress, 
OLMax is the normalized Maxwell stress, and the double 
angle brackets (()) denotes both time and space averages. 
However, there is no guarantee that a is a constant in 
disks with large scale density structures (e.g. gaps). Nor- 
mally, global simulations (Sorathia et al. 2012) find large 
a fluctuations at different parts of the disk. 

Gap opening by massive planets (much larger than a 
thermal mass) in MRI-turbulent disks has been stud- 
ied by several authors using net toroidal flux or zero 
flux MHD simulations (Winters et al. 2003, Nelson & 
Paploizou2003, Papaloizou et al. 2004, Uribe et al. 2011). 
Global simulations by Nelson & Papaloizou (2003) sug- 
gest the gaps in MRI disks are slightly deeper and wider 
compared with those in viscous disks, while Winters et 
al. (2003) found that gaps are slightly shallower in MRI 
disks. Overall, they found that the shapes of gaps opened 
by high mass planets are not very different between MRI 
and viscous disks. 

Papaloizou et al. (2004) found that both local shear- 
ing box MHD simulations and global MHD simulations 
show good agreement on the properties of gaps. The 
local approximation enables high resolution simulations 
to bo carried out. Furthermore, local shearing box sim- 
ulations can maintain a steady state, which is crucial 
to study planetary wake and gap properties since some 
features can only be identified by averaging over hun- 
dreds of orbits. Although local simulations lack curva- 
ture terms which prohibit studying the differential torque 
which controls planet migration, it is still possible to 
study the one-sided Lindblad torque and the gap open- 
ing process in the local shearing box approximation as 
long as the gap structure is sharper than the disk radius 
R so that the curvature terms are not important. 



The tidal torque felt by the planet can be affected by 
the presence of magnetic fields in disks, which has been 
studied by linear calculations (Terquem 2003) and nu- 
merical simulations (Fromang et al. 2005, Muto et al. 
2008). However, these calculations assume strong mag- 
netic fields and non-turbulent disks. It is unclear how 
MRI turbulence will affect the tidal torque. Recent 
global MHD simulations have shown some additional 
torque close to the planet (Uribe et al. 2011; BarutcaTi ct 
al. 2011). Local shearing box limit is a good way to test 
if this additional torque is associated with disk global 
curvatures. 

Most previous work on gap opening in MHD disks as- 
sumed zero net magnetic flux in simulations. In this 
paper, we focus on MHD disks with a small net verti- 
cal magnetic field. The reason for choosing a small net 
field is three-fold. First, recent unstratified simulations 
have shown that the turbulent properties depend on both 
the numerical resolution and box size with zero net field 
(Fromang & Papaloizou 2007). Second, simulations with 
net magnetic flux are more realistic. Wc expect that 
small local patches of the disk will be threaded by net 
vertical fields from molecular clouds, or by toroidal fields 
generated by the large scale disk dynamo, or even both 
(Sorathia et al. 2010). Third, by varying the net field 
strength, we can control the strength of the MRI tm-bu- 
lence, which allows us to systematically study gap open- 
ing in disks with various turbulence stresses. 

Furthermore, previous work which study gap open- 
ing focused on deep gaps induced by massive planets 
("Jupiter" mass range). Here we are interested in cases 
where gaps just start to open in order to test the "ther- 
mal" and "viscous criterion". This motivates us to fo- 
cus on gap opening by low mass planets ("earth" mass 
range). 

This paper is organized as follows. In §2, we intro- 
duce our methods. To provide a benchmark for compari- 
son, we first study gap opening in inviscid hydrodynamic 
(HD) disks in §3, while planetary wake properties and 
gap opening in MRI-driven turbulent disks is presented 
in §4. Gap opening criteria are discussed in §5. After 
discussions in §6, we conclude the paper in §7. 



2. METHOD 

2.1. Basic Equations 

This study uses Athena (Stone et al. 2008), a higher- 
order Godunov scheme for MHD with piecewise parabolic 
method (PPM) for spatial reconstruction, the corner 
transport upwind (CTU) method for multidimensional 
integration, and the constrained transport (CT) to con- 
serve the divergence-free property for magnetic fields 
(Gardiner & Stone 2005, 2008). The simulations are set 
up as isothermal unstratified disks in the local shearing 
box approximation (Hawley et al. 1995, Stone & Gar- 
diner 2010), in which the MHD equations are solved in 
a reference frame centered at radius i?o corotating with 
the disk at orbital frequency fio = ^{Rq)- Ignoring cur- 
vature terms, the isothermal MHD equations are 



|+V.(pv)=0, 



(4) 
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^ + V-(pvv + TB + Tv) + Vp = 2qxpr2^i-217okxpv, 
ot 

(5) 
(6) 



— - Vx(vxB) = 



where p is the gas pressure, the magnetic stress tensor 
Tb is 

Tb = (BV2)I-BB, (7) 



and the shear parameter q is 

1 dlnOF- 

q 



2 dlnr 



(8) 



so that for Keplerian flow q=3/2. Note that these equa- 
tions are written in units which assume the magnetic 
permeabihty /i =1. An equation of state for an isother- 
mal gas p = c^p is used. 

Note that, for HD disks, we have included explicit vis- 
cosity in the momentum equation, through the viscous 
stress tensor Tv 



(9) 



Tv;ij = I'P [ diVj + djVi - -5kVk(5ij 



where v is the kinematic viscosity. This allows us to 
study HD models with viscosity, and compare these to 
the MHD models which are always inviscid. 

An orbital advection scheme (Masset 2000) has been 
implemented in Athena (Stone & Gardiner 2010). The 
y-direction momentum equation is split into a linear op- 
erator for advection with the background flow velocity, 
as well as the MHD equations for the velocity fluctua- 
tions. The orbital advection scheme significantly accel- 
erates the simulation since the maximum allowed time 
step is limited by the magnitude of the velocity fiuctua- 
tions instead of the orbital velocity. It also reduces the 
inhomogeneous truncation error produced by differential 
rotation (Johnson et al. 2008). 

Gressel & Ziegler (2007) have pointed out that the 
shearing box boundary conditions can destroy conser- 
vation, since the integral of the fluxes over the two radial 
faces are not identical to machine precision. However, in 
order to model the density structure (gap formation) in 
the disk, it is important that the total mass is conserved 
to round-off error. Thus, we have used the remap scheme 
similar to Stone & Gardiner (2010) for the density fluxes. 
First, we remap the radial component of the mass flux at 
each radial face. Then, we use the arithmetic average of 
the radial mass flux computed for each grid cell at each 
radial face, and the remapped value of the radial mass 
flux from the corresponding grid cell on the opposite face, 
to update the density in the cells next to the boundary. 
This scheme guarantees the integral of the mass flux is 
identical at two boundaries in the x-direction. 

To simulate the effect of a planet in 3D unstratifled 
shearing box simulations, a line mass (symmetric in the 
z— direction) is placed at the center of the box. This 
geometry is identical to the symmetry of 2D HD sim- 
ulations, but allows MHD turbulence to be modeled in 
full 3D. For comparison to 2D simulations we average 
quantities in the z— direction. 

The planet potential is smoothed over several grid cells 
to avoid small time steps associated with the divergence 



of a point source. Following Dong et al. 2011a, we use 
fourth order smoothing 



-GM, 



ml/2 



(10) 



where d is the distance to the z-axis [d — \/ + y"^) and 
Rs is the smoothing length. This potential converges to 
that of a point mass as {Rg/dY for d^ Rg and deviates 
from the point mass potential by only 1% at 2.3 Rg. 

Due to the periodic boundary condition, the gradi- 
ent of the potential has a discontinuity at the boundary. 
Thus we smooth the potential as 



$s,p(rf) = %{df) - ^{%{df) - %{d)Y + G^MIIRI, 
\id<df (11) 

and 

$,,p(d) = $p(d/) - GM^jRgo iid> df (12) 

where d/ is the cut-off radius, beyond which the potential 
is flat, and Rgo is a smoothing length at the cut-off radius. 
For large Rso, eq. (dJ) reduces to pUj) . With a flnite 
Rso the potential becomes flat smoothly (the gradient is 
a continuous function) at o? = -Rso- In all our simulations 
df=7.5 H and Rso=5Q H, where H is the disk scale height. 
With this set-up this potential differs from a point mass 
potential only by 3.5% at df—6 H. 

Introduction of the planet potential can disturb the 
background flow in the disk signiflcantly, and with pe- 
riodic boundary conditions this initial disturbance can 
affect later disk evolution. Thus, in HD simulations, we 
linearly ramp up the amplitude of the planet potential 
from to 5 orbits, while in MHD simulations, we linearly 
ramp up the potential from 5 to 10 orbits after the MHD 
turbulence is established at 5 orbits. 

The smoothing length in eq. (|10p for simulations with 
a 0.1, 0.3, and 1 Mj planet is chosen as 0.125, 0.18, and 
0.269 H respectively, so that the free-fall timescale onto 
the planet (dcflned as Eq. 29 of Dong et al. 2011a) is 
0.125 in our code units. The numerical time step is 
flxed at 2x10^"^ ('^60 times smaller than the free fall 
timescale) even though the orbital advection scheme al- 
lows for a much longer time step. This small time step 
allows us to trace the fluid motion around the planet ac- 
curately enough to avoid unsteady fluid motion (Dong 
et al. 2011a). With Athena, Dong et al. (2011a) found, 
for the case they studied (with resolution 64 grid points 
per scale height), that the numerical time step needs to 
be at least 140 times smaller than the free fall timescale 
onto the planet. While it is clear from physical grounds 
the time step must be smaller than the free fall time near 
the planet to properly resolve orbital motion of gas near 
to the planet, the exact ratio between the time step and 
the free fall time can vary between different codes (Kley 
et al. 2012). By exploring a larger parameter space, we 
have also found that the critical time step is roughly pro- 
portional to the grid size, as one might expect from the 
CFL condition. Thus, with a factor of 2 lower resolution 
in our simulations (32 grid points per scale height as 
discussed below), the critical time step can be 70 times 
smaller than the free fall timescale. We have run HD 
simulations with our choice of parameters and conflrmed 
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that all the runs are stable for 400 orbits with the same 
boundary set-up in Dong et al. (2011a). 

2.2. Model Set-up 

We have carried out both 2-D HD (inviscid and vis- 
cous) and 3-D MHD simulations with a 0.1, 0.3, and 1 
thermal mass {Mth) planet at the box center. The code 
units are chosen as Cs = 1, $7 = 1, and po = Ij so that 
the disk scale height (H) and the thermal mass {Mth) 
are both 1. The numerical resolution is 32 grid points 
per H, and the box size is 16 H xl6 H x 1 H in most 
simulations. The parameters of the simulations are sum- 
marized in Tables 1 and 2. We divide the simulations 
into one of the three categories. 

1) Inviscid 2-D HD simulations with various box sizes 
and a 0.1 Mt^ planet in the disk. This set of simula- 
tions is designed to test whether linear waves launched 
by the planet will steepen into shocks within the sim- 
ulation domain, and whether shocks entering from the 
radially periodic boundary conditions significantly affect 
the gap region. We use these simulations to determine 
the ideal box size for further study. 

2) 3-D unstratified ideal MHD simulations with and 
without planets. Simulations without planets allow us 
to study density fluctuations and turbulent stress due 
only to the MRI. These calculations serve as a reference 
for comparison to simulations with planets in order to 
distinguish, for example, the density dip due to a gap 
opened by a planet from that associated with zonal flows 
due to the turbulence (Johansen et al. 2009, Simon et al. 
2012). 

MHD simulations with different mass planets (0.1, 0.3 
and 1 Mf) have been studied. Two different initial 
vertical field strengths have been applied with /3o=400 
and 1600 (the ratio between the initial gas pressure and 
the magnetic pressure), leading to disk MHD turbulence 
with a= 0.17 and 0.085. The initial net vertical flux 
is conserved in these simulations. In ideal MHD, the 
wavelength for the fastest growing linear MRI mode is 

H/\ = 0.109/3^^ (Hawley et al. 1995). With /3o=400, 
and 1600, our vertical box size (H) fits 2 and 4 wave- 
lengths of the fastest growing mode. With 32 grids per 
H, the most unstable wavelength is resolved by 16 and 
8 grid cells respectively. This is sufficient to ensure the 
properties of the MRI-driven turbulence are numerically 
converged (Hawley et al. 2011, Simon et al. 2009, So- 
rathia et al. 2012). 

3) 2-D viscous HD simulations with constant kinematic 
viscosity that have the same R — (j) stresses as those in 
MHD simulations. The gap structures in these viscous 
disks are compared with the gaps in the MHD simula- 
tions from 2) above. The planet masses are again 0.1, 
0.3 and 1 Mth- Note that in viscous disks the xy com- 
ponent of the viscous stress tensor Ty ~ v p{dvy / dx + 
dvx/dy) ~ lypol.^il. Equating r„ with the MRI stress 
Tmb.1 ~ apoCj leads to v — acl/l.bft. For an exam- 
ple, in order to compare with the MHD run M10B400 
which has a=0.17, the viscous run VIO needs to have 
zy=0.17/1.5=0.113. 

2.3. Boundary Conditions 

In our simulations, we have chosen the periodic bound- 
ary condition in both the x— (radial) and y— (azimuthal) 



directions. A periodic boundary condition is natural in 
the azimuthal direction (y or 9). However, applying it to 
the radial direction effectively introduces an infinite grid 
of virtual planets at integer spacings of the radial box 
size. A more natural radial boundary condition would be 
outflow, or the wave damping condition used by de Val- 
Borro et al. (2006). However, to model MHD turbulence 
in the local shearing box periodic boundary conditions 
must be used at both x and y directions. Therefore, 
in order to study the effect of the x-direction periodic 
boundary condition on our results, and to choose the 
optimum domain size, we have calculated three inviscid 
HD disks (runs II, 12, 13) having different radial extents: 
8Hxl6H, 16Hxl6H, and 32Hxl6H and a 0.1 Mf planet 
at the center. These simulations will be discussed in de- 
tail in §3, and in the following paragraph we focus on the 
effect of the boundary condition. 

Figure [T] shows the disk surface density contours at 
400 orbits for all three simulations (runs II, 12, 13). The 
effect of the x-direction periodic boundary condition is 
most apparent in the smallest box (8 H x 16 H). The 
wakes from nearby "virtual" planets overlap with the 
gap region in the box. Multiple stripes appear in this 
box instead of two clean gaps, and the edge of the gap is 
close to the a;-direction boundaries. While, the biggest 
box (32 Hx 16 H) has the cleanest gap, the intermediate- 
sized domain (16Hx 16 H) resolves the two gap struc- 
ture well, with the two gaps well within the simulation 
domain. Furthermore, the wakes from nearby planets 
are much weaker than the wake from the planet at the 
center. Thus, as the best compromise between cost and 
accuracy, we have chosen 16H x 16Has our box size 
for most other simulations. We note that in both vis- 
cous and MHD models presented later, the wakes from 
"virtual" planets are a lot weaker than those in the invis- 
cid HD model presented here since the wakes are quickly 
damped either by viscous stress or the MRI turbulence 
(§4). The box size used in this study is therefore even 
more suitable in these cases. 

3. GAP OPENING IN INVISCID HD DISKS 

Figure [1] also questions the "thermal criterion" for gap 
opening in the inviscid fluid limit. Gaps are opened by a 
planet with mass significantly smaller than the "thermal 
mass" . The objection to the "thermal criterion" was first 
suggested by Goodman & Rafikov (2001), stating that 
linear waves from a low mass planet (Mp < Mth) steepen 
to shocks after they propagate for a distance 



\Xsh 



0.93 



7 + 1 Mp 
12/5 Mth 



-2/5 



H . 



(13) 



where 7 is the adiabatic index. Shocks efficiently transfer 
the wave angular momentum to the background flow, and 
two gaps should gradually open at the distance of \xsh\ 
on each side of the planet. Thus the thermal mass Mth 
is not a gap opening threshold in this nonlinear wave 
propagation theory. For a Mp = Q.lMth planet in an 
isothermal disk (7=1), the waves shock at the distance 
-2.5ff. 

In Figure [11 we clearly see gap opening for a 0.1 Mth 
planet in an inviscid HD disk (runs II, 12, 13 ). The gaps 
are indeed deepest at \x\ Xsh 277. The gaps are more 
prominent with bigger boxes since they are less affected 
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by the boundary conditions, as discussed in §2.3. 

Figure [2] shows the space time plot for the disk surface 
density in run 12 (the same simulation as in the middle 
panel of Figure [ij . The x-axis is the time in units of 
the orbital period after the planet is inserted while the 
y-axis is the y-direction (azimuthal) averaged disk sur- 
face density across the gap. We can see that two gaps 
are gradually opened starting at ~ 2i?, getting deeper 
and wider with time. For a long time (hundreds of orbits) 
a ribbon of material is maintained between the gaps, in 
agreement with the predictions of Rafikov (2002b) . This 
structure appears because density waves do not dissipate 
in the coorbital region in an inviscid case, and thus an- 
gular momentum is not transferred to the disk, leaving 
the fluid between the gaps essentially unperturbed. The 
symmetry of the two-gaps structure in our simulations is 
due to the shearing box geometry, and should be broken 
in global simulations, as well as if planet is migrating 
through the disk. 

No steady state is reached, since in an inviscid HD disk, 
there is no torque to balance the torque from shock dissi- 
pation. We find that gaps continue to widen and deepen 
until some instability, possibly the Rossby wave instabil- 
ity (Lovelace & Hohlfcld 1978; Lovelace et al. 1999; Li 
et al. 2000; de Val-Borro 2007; Lin & Papaloizou 2010), 
or Rayleigh instability sets in and destroys the ribbon of 
gas coorbital with the planet. This should finally result 
in a single gap in a disk, but reaching this state can take 
long time for a low mass planet in an inviscid disk. 

The timescale to open a factor of 2 deep gap is ~400 
orbits in our local simulations. However, to derive a 
timescale for gap opening in a realistic disk based on 
our local shearing box simulations, the time in our sim- 
ulations needs to be multiplied by 2ttR/16H , where R 
is the distance from the central star to the planet and 
16H is our box size in y-direction. This is simply due 
to the fact that the angular momentum injected by the 
planet to remove all the material per unit radial distance 
is proportional to the box extent in the y-direction (which 
should be 2ttR instead of 16H), while the torque ex- 
cited by the planet and the amount of shock dissipation 
per unit radial distance are independent on the box sizes 
0. Thus, for a global disk with H/R=0.1, the timescale 
to open a gap by a 0.1 Mth planet is estimated to be 
400x27ri?/16i?- 1,500 orbits. 

4. GAP OPENING IN MRI DISKS WITH NET VERTICAL 

FLUX 

In this section, we will first study the wake properties 
and wake damping mechanism in MHD turbulent disks 
using simulations with low mass planets (e.g. M01B400 
and M03B400, §4.1). Then we will compare the gap 
properties between the MHD turbulent disks and the vis- 
cous disks using simulations with more massive planets 
(e.g. M10B400, M10B1600, §4.2). 

The disk surface densities from these simulations are 
shown in Figure [3] for planet masses of 1 Mf, 0.3 Mt, and 
0.1 Mt (from top to bottom). The left column shows 
a snapshot in the MHD simulations (runs M10B400, 
M03B400, and M01B400) at 70 orbits. The turbulent 
nature of the MRI disk can be clearly seen. Density 

^ The box needs to be larger than ~ 4H X AH so that most of the 
planetary torque has been excited (Goldreich & Tremaine 1979). 



waves excited by the MRI turbulence (Heinemann & Pa- 
paloizou 2009) are present across the whole simulation 
domain. The standard deviation for density fluctuations 
relative to the background density due to the MRI tur- 
bulence alone is as large as ■^O.S. With such large density 
fluctuations, both the gap structure and planetary wake 
are hard to discern. In order to extract the gap and plan- 
etary wake from MRI turbulence, we have averaged the 
disk surface density every orbit from 80 orbits until the 
end of the calculation (shown in Table 2). These aver- 
aged surface densities are shown in the middle column 
of panels, where both the gap structure and planetary 
wakes are clearly evident. By comparing to gaps in vis- 
cous HD disks (runs VIO, V03, VOl, the right column of 
panels), it is clear that the gaps in MRI disks are deeper, 
which will be further explored in §4.2. 

It is important to note that MRI-turbulent disks can 
produce large scale density structures that persist over 
tens to hundreds of orbits, called zonal flows (Johansen 
et al. 2009, Simon et al. 2012), even without the pres- 
ence of a planet. To distinguish the gap opened by the 
planet from zonal flow, we show the space-time plot for 
the y-averaged disk surface density for runs M10B400, 
M03B400, M01B400 and B400 in Figure IH The magni- 
tude of density fluctuations produced by zonal flow can 
be seen in the bottom panel for the MRI disk with no 
planet case(B400). Clearly, the density dip in the 1 Mt 
planet run (M10B400, the top panel) is much larger in 
magnitude than the zonal flow, and the position of the 
gap is close to the planet. Thus we can confidently say 
that it is indeed a gap induced by the planet in the 
M10B400 case. Before we study the gap opening pro- 
cess in MHD disks (§4.2), we first need to understand 
how waves deposit angular momentum in MRI turbulent 
disks (§4.1). 

4.1. Wake Properties 
4.1.1. Identifying Wakes 

The damping mechanism of propagating waves deter- 
mines how angular momentum is deposited into the disk 
and furthermore the gap opening process. Various damp- 
ing mechanisms have been proposed, including viscous 
damping (Takeuchi et al. 1996) in viscous disks and 
shock damping (Goodman & Rafikov 2001) in inviscid 
HD disks. In this work we will investigate the damping 
mechanism in the high Reynolds number MHD turbu- 
lence expected in real protoplanetary disks. 

To identify the shock region, in Figure [5] we have 
plotted the vertical average of the velocity divergence 
r/ = (Aa;/cs)V • v (using the sound speed and grid spac- 
ing to normalize the velocity and length scales respec- 
tively) for both MRI disk (B400, bottom panels) and 
MRI-fplanet disk (M10B400, top panels). In our simula- 
tions, we find 77 < —0.2 is a good criterion for identifying 
strong shocks. 

The leftmost panels of Figure [5] show the contours of 
77 at 80 orbits. The next column of panels outline the 
region where 77 <-0.2, while the next column after that 
show these outlines in every snapshots from 80 to 130 
orbits using a 5 orbit interval. Even in the case of no 
planets (B400), as shown in the bottom panels, sharp 
density features are present in the MRI-driven turbulent 
disks. These shocks in MRI disks are consistent with 
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structures explored in Heinemann &; Papaloizou (2012). 
However, due to the random nature of the turbulence, 
these shocks appear uniformly in the disk (the lower 
right panel). When a 1 Mth planet is present in the 
disk (M10B400), as shown in the upper panels of Figure 
[51 it excites strong shocks around the planet, and rj can 
be as low as -0.7. For comparison, a planet in a 2-D in- 
viscid HD disk (IlOinv) also excites a strong shock (as 
shown in the upper rightmost panel), and the smallest 
rj is also -0.7. On the other hand, in the equivalent vis- 
cous disk simulation (VIO), rj of the wake is only -0.066, 
one order of magnitude smaller than that in both the 
MHD and inviscid HD cases. Thus, we conclude that 
shock damping is the main wake dissipation mechanism 
in MRI disks, which is significantly different from viscous 
damping with a kinematic viscosity. 

As shown in the upper panels of Figure [SJ the wake's 
shock position in MRI-turbulent disks varies with time 
due to the disk turbulence (the second to right panel), 
which is different from the static shock in inviscid HD 
case (the rightmost panel). The random fluctuations in 
the shock position appear to spread it out over a finite 
distance, even though at any instant in time the shock 
profile is sharp. As we will show below, the time aver- 
age density profile of the planetary wake is surprisingly 
similar to the wake in a viscous disk. 

4.1.2. Averaged Wake Density Profile 
The one-armed planetary wake, which is located at 

3a;2 
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in the local limit (Goodman & Rafikov 2001), is due to 
the coherent interference (Ogilvie & Lubow 2002) of all 
the modes excited by the planet (perturbed quantities 
oc exp[i{m(j) — ujt) + kuR)]), each satisfying the WKB 
dispersion relationship without self-gravity 
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where k is the epicyclic frequency {k = n for a Keplerian 
disk), kn is the radial wavenumber in the WKB limit, 
and rip — n{Rp) is the pattern speed io/m in the inertial 
frame. 

Since this one-armed wake is the "signpost" of planet- 
disk interaction, it is worth studying its structure in de- 
tail. The density profiles across the wake in the invis- 
cid HD disk have been studied by Goodman & Rafikov 
(2001), Dong et al. (2011a) and Rafikov & Petrovich 
(2012). In this work, we will examine the density profiles 
of the wake in MRI-driven turbulent disks. In a single 
snapshot, the planetary wake is significantly disturbed 
by MHD turbulence (as shown in §4.1.1 above). Thus 
we average the disk surface density at every time step 
over 100 orbits for the planetary wake to stand out. 

In order to compare the wake at different x positions 
in different simulations, we need to apply some scaling 
procedures as in Goodman & Rafikov (2001). For HD 
simulations, the wake density profiles are measured along 
azimuthal cuts through the disk at a:=0.5, 1, 2, 4 H, and 
then we shift the derived one dimensional density profiles 
in y by yt given in Eq. We then scale the den- 

sity to a dimensionless quantity as follows (Goodman & 
Rafikov 2001): first, the mean azimuthal density at x is 



subtracted to derive JE. Since increases with x in the 
linear regime due to the angular momentum flux conser- 
vation, we normalize ST, by x^^"^. Finally we scale ST by 
T{x){Mp/Mth), where S(x) is the y-direction averaged 
surface density along x. In the linear theory, with this 
scaling, the density profile should be independent of Mp 
at a fixed separation x, and the magnitude of the wake 
density profile should remain constant with increasing |a;| 
(until the wake shocks). 

For MHD simulations, the above scaling procedure is 
similar but the scaling speed now is the fast magnetosonic 
speed y^cl + v\, where va is the Alfven speed, instead of 
the sound speed {cg) as in HD cases (Terquem 2003). We 
find that the MRI saturated states in these simulations 
have (i «4, thus the fast magnetosonic speed is 1.22 c^ 
and the new length scale is 1.22 H. Therefore we make 
the density cuts at x=(±)0.61, 1.22, 2.44, 4.88 H. Then 
we shift the one dimensional density profile by 
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similar to Papaloizou et al. (2004). Finally we scale ST by 
To{Mp/M[^) where M^'^ = (c^ -|- v\f''^l{Gn) is a new 
thermal mass defined via the fast magnetosonic speed. 
Since the resulting density profiles have a slight asymme- 
try due to zonal flow in the disk, we average the density 
proflles from both the positive and negative x sides of 
the disk. 

The resulting wake density profiles from a 0.3 Mt 
planet in disks are shown in Figure [51 Solid, dotted and 
dashed curves are from MHD simulations (M03B400), 
viscous HD (V03), and inviscid HD simulations (I03inv) 
respectively. Remarkably, with the new scaling Eq. (jl6p , 
the averaged density profiles from the MHD case are 
quite similar to the density profiles in the viscous disk 
with the same stress. The sharp density profiles in the 
inviscid HD disk (I03inv) indicate the formation of the 
strong shock, while the smooth profiles in the viscous 
case (V03) indicate the viscous damping. However, the 
smooth density profiles in the MHD case suggest that, 
although the damping mechanism is the shock damping 
(§4.1), the temporal fluctuation of the shock position in 
an MRI-turbulent disk spreads the wake and, in a time- 
averaging point of view, the shocks in turbulence have 
similar effects as the viscosity on the averaged wake den- 
sity profiles. 

4.1.3. Averaged Torque 

Using the time averaged surface densities (E) from 
MHD simulations, we can calculate the time averaged 
torque density (dT/dx) from the planet by integrating 
{T)d(t)/dy over the y-direction 



dT f/^.d(j). 



(17) 



where the cylindrical region centered on the planet within 
the planet's Hill radius is excluded for the torque calcu- 
lation since this region should be bound to the planet 
and co-moving with the planet. Note that we can cal- 
culate the torque density by only integrating (E) over 
y direction instead of integrating (p) over both y and z 
directions, since the planetary potential is independent 
of the z-direction in our calculations. 
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The averaged torque densities for M10B400, M03B400, 
and M01B400 are shown in Figure [7] as the sohd curves. 
By comparison, the dashed curves are from inviscid HD 
simulations (IlOinv, lOSinv, lOlinv), while the dotted 
curves are from viscous HD simulations (VIO, V03, VOX). 
As well known, the torque densities for inviscid HD disks 
(dashed curves) have peaks around 1 H, due to the torque 
cutoff (Goldreich & Tremaine 1980, Artymowicz 1993, 
Rafikov & Petrovich 2012). In contrast to both the invis- 
cid and viscous HD cases, the averaged torque densities 
from MRI disks have larger peak values at smaller |a;|. 
More discussions are in §7.1. 

4.2. Deeper and Wider Gaps 

As shown in Figure[31 the gaps in MHD-turbulent disks 
with net vertical flux (M10B400) are significantly deeper 
and wider than the gaps in viscous disks with the same 
stress (VIO). To be more quantitative, we average the 
disk surface densities from MHD simulations M10B400 & 
M10B1600 over both time, and the y— and z— directions 
to derive the surface density profile along x, and compare 
them with those from viscous simulations VIO & VlOsv, 
as shown in Figure |S1 Due to the large mass concentra- 
tion in the circumplanetary region around the planet, the 
region within \y\ < H is removed for the averaging. The 
solid curves are for MHD simulations, while the dotted 
curves are for viscous HD simulations. For run M10B400 
(solid black curves), the ratio between the smallest sur- 
face density within the gap and the surface density at 
the box edge is ~ 50%, while this ratio is ~ 80% for 
the viscous run VIO (dotted black curves). With smaller 
stresses, the gaps in M10B1600 (solid grey curves) and 
VlOsv (dotted grey curves) are both deep, but the gap 
in the MHD case M10B1600 is still significantly deeper 
and wider than the corresponding viscous case VlOsv. 

In order to understand why the gaps are significantly 
deeper and wider in these MHD cases, we note that the 
steady-state gap shape is determined by the balance be- 
tween angular momentum deposited into the disk and 
the stress from MRI turbulence (or viscosity). For vis- 
cous HD disks, this balance can be expressed as 

^"/"j J d J {pv^Vy)dydz d J -upid^Vy + dyV^)dydz 

p-^dydz —2 = 5 2 

ay dx dx 



while for the MHD disk it gives 
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while for MHD disks, it is 
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Note that we separate the Reynolds stress in Eq. (|19p 
into two parts: the portion excited by the planet and 
that excited by MRI turbulence. 

Next, we integrate the above equations along x start- 
ing from the planet position x=0, and then perform an 
ensemble average for both time and the z-direction. For 
the viscous HD disk this gives 

i-^)'^^ - / {{pvxvy))dy 



{{vpidxVy + dyVx)))dy 



(20) 



For the disk region from to x, the left-hand sides of Eq. 
(120]) and ([2T|) represent this region's net gain of angular 
momentum from the waves excited by the planet. They 
consist of the integrated torque exerted by the planet 
(the first term), minus the angular momentum fiux of 
the waves leaving this region (the second term). The 
right-hand sides of Eq. pU)) and (|2T|) represent the net 
flux of angular momentum into this region carried by 
viscosity, or by MHD stresses. A steady gap shape is 
achieved when the two sides of each equation balance. 

We integrated the torque density from the same mass 
planet for both MHD and viscous runs (as derived in 
§4.1.3), and found that the integrated torques (the first 
terms on the left-hand singes of Eqs. [20land [2l]) are close 
to each other within 20% difference. Considering the 
torque density normally is much larger than the planet 
induced Reynolds stress, the left hand sides of Eqs. (|20p 
and ((2T|) are similar. This suggests that the rate of depo- 
sition of angular momentum from the planet to the disk 
is similar in both MHD and viscous HD disk simulations, 
and therefore the differences in gap depth between these 
two cases can only be due to the stress gradient along 
X (the right hand sides of Eqs (|20p and (I2T]) ') as shown 
below. 

4.2.1. Viscous versus Maxwell stress 

The averaged viscous stresses for the viscous disks 
(VIO, VlOsv) and the Maxwefl stresses for the MHD 
disks (M10B400, M10B1600) are shown in the upper 
right panel of Figure [5] Only xy components of the 
stresses are shown. By comparing with the density pro- 
files in the upper left panel, we see that the viscous 
stresses are almost proportional to the densities in the 
viscous disks (dotted curves). This correlation is ex- 
pected since the shear in the y direction dominates the 
stress and the stress can be approximated by l.bi'pQ. v 
and O are both constant in viscous disks. However, as 
shown as the solid curves, in MHD disks, the Maxwell 
stresses across the gaps are significantly smoother than 
the density changes. Since the gap shape is not deter- 
mined by the stress itself but by its gradient along x (Eq. 
[TO]) , this smoother stress with respect to p across the gap 
means that to maintain the same stress gradient to bal- 
ance the planetary torque the gap needs to be deeper in 
the MHD case. 

Along the x direction, the smoother profiles of the 
Maxwell stress compared with the density can also be 
characterized by dividing the stress with density, or 
equivalently calculating a parameter (Eq. [3]). For the 
viscous cases, a is almost constant across the gap, as 
shown in the upper right panels of Figure |9l and the 
time and yz-direction averaged a is shown as the dot- 
ted curves in the lower left panel of Figure [S) For MHD 
disks with gap opening, (Xuax increases towards the gap 
region, as shown in the lower right panel of Figure [9] and 
the averaged auax is shown as the solid curves in the 
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lower left panel of Figure [SJ again suggesting that the 
MHD stress is smoother across the gap than the density 
or viscous stress. The time and yz-direction averaged 
magnetic field (net Bz) is also stronger within the gap 
(the lower right panel of Figure [8]) . 

We also explored cases with different planet masses 
(e.g. 3 Mt, 0.3 Mt), and box sizes (e.g. 32 H in y di- 
rection) , as shown in Figure [TOl aMax in these cases 
all show the same trend of increase towards the gap re- 
gion, and gaps are again deeper and wider in MRI disks 
than corresponding viscous cases. The density and a 
structure in the xy plane for our biggest box simulation 
(M10B1600b) are shown in Figure [TI] which nicely illus- 
trate the turbulent nature of the disk and demonstrate 
the higher a in the gap region. 

In order to quantify this weak stress dependence on the 
density across the gap, We calculate the deviation of the 
Maxwell stress along x 

{{BxBy))a {{BxBy))o 

with respect to the deviation of the density 
(A((/9))/((p))o), as shown in the left panel of Figure 
[ll Both {{B^By)) and {{p)) are derived from Figure H 
and Figure [To] and {{BxBy))o corresponds to {{p))o = 1. 
For the viscous stress (T^ — l.bvpVl), ATy/T^j = Ap/p. 
However, for all the MHD cases with net vertical 
magnetic flux, A{{B,By)) / {{B,By))o = CA((p))/((p))o 
where C - 0.25 < 1, suggesting {{B^By}) oc {{p))°-^^ 
across the gap, and indicating the weak dependence of 
the stress on the density across the gap. Thus, it is 
possible to use this new relationship between stress and 
density to construct a viscous disk model to simulate 
the gap as in the MRI cases. However, we caution that 
this relationship is only approximately true for these 
limited cases as shown. 

Finally, in order to understand why aMax is larger in- 
side the gap, we test if the empirical relationship between 
the turbulent Maxwell stress aMax and the plasma ((/3)) 
in MRI saturated states {aMax ~ 1/2 ((/?))) (Hawley et 
al. 1995) still stands across the gap locally. This is shown 
in the right panel of Figure 1121 where each point repre- 
sents a y direction averaged aMax with respect to the 
averaged at each radius x. As x marches towards 

the center of the gap, the aMax value increases and the 
{{(3)) value decreases. And they tightly following the 
empirical relationship. The decreasing towards the 
gap suggests that the turbulent magnetic fields are more 
concentrated in the gap region with respect to the den- 
sity. 

Not only the turbulent magnetic fields, but also the net 
magnetic fields concentrate to the gap region (the lower 
right panel of Figure (S]) . These two evidences suggest 
efficient global transport of magnetic fields to the low 
density gap region in MHD disks with net vertical fields. 

4.3. Persistence of the gap 

As discussed in §4.2.1, the effective a is higher in the 
gap region in net vertical flux MHD disks. Thus, if the 
gap depth is similar between MHD and viscous disks, the 
Maxwell stress in net vertical flux MHD disks will have 
a smaller gradient than the viscous stress in viscous HD 
disks. Since it is the stress gradient that determines the 



global transport of the angular momentum and mass, any 
pre-existing density features will behave quite differently 
among MHD disks and viscous HD disks. 

To test this hypothesis, we restart runs M10B1600 
and VlOl at a time of 50 orbits, when gaps formed by 
the planet are fully developed, but with planetary po- 
tential switched off in the subsequent evolution. The 
space-time plot for the disk surface densities are shown 
in Figure [131 As clearly shown in the upper panel, the 
gap feature in M10B1600 even persists for 50 orbits after 
the planet disappears, while in viscous disks the gap is 
quickly closed for ~10 orbits comparable to the viscous 
timescale {2HY/v ^ 20 orbits. This confirms that the 
stress is uniform in net vertical flux MRI disks, leading to 
the inefficiency in erasing the preexisting density features 
in the shearing box. 

This result has far reaching implications for the mea- 
surement of the "viscous" timescale from the diffusion of 
density features in observed disks. If the disk is threaded 
by a vertical magnetic field, the disk evolution timescale 
will be longer than the "viscous" timescale. 

5. GAP OPENING CRITERIA 

5.1. "Thermal Criterion" 

We have demonstrated that in an inviscid HD disk, 
a planet with a mass 10 times smaller than the ther- 
mal mass can open gaps by shock formation. Unlike gap 
opening by a massive planet (> IMtu) where a single gap 
is opened around the planet, a low mass planet opens two 
gaps away from the planet at the position where density 
waves shock (Eq. [T3l) . 

Simulations with a 0.3 Mth planet embedded 
(M03B400, M03B1600) start to show gap features around 
the planet (lower panels in Figure [TU)) . Similar to 1 Mth 
mass cases, a smaller aMax leads to a deeper gap. We 
expect a 0.3 Mth planet can open a clearer gap in a even 
less turbulent disk. Thus, we do not confirm the sugges- 
tion of Papaloizou et al. (2004) that the thermal mass is 
a threshold for gap opening. Low mass planets can also 
open gaps as long as the disk stress is low. 

On the other hand, the timescale for gap opening is 
longer for smaller planets, and eventually the feedback 
due to migration starts playing an important role in the 
gap opening process(Ward 1997, Rafikov 2002), which 
requires global simulations (Li et al. 2009). 

In a future paper, we will show that there is also clear 
evidence that the gap structure not only depends on the 
absolute value of the stress (a) but also on the field ge- 
ometry. 

5.2. "Viscous Criterion" 

The derivation of the traditional "viscous criterion" 
(Eq. [2]) assumes the thermal criterion is satisfied so 
that the disk scale height is replaced by the planet's 
Hill radius. If we relax this assumption, and balance 
the momentum flux excited by the planet (Goldreich & 
Trcmainc 1980) 

Fff «0.93(GMp)2 5^^^ (23) 

and the momentum flux due to viscosity 

= iirY^vRln (24) 
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the viscous criterion becomes 
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where 6 is just an order of magnitude estimate and the 
detailed value depends on how we define a "gap" . How- 
ever, this criterion needs to be modified to study gap 
opening in the shearing box limit, since the momentum 
flux due to viscosity in the shearing box is 



(26) 



where Hy is the box size in y direction. Thus, the viscous 
criterion in the shearing box is 
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For a disk having i/=0.113, the viscous criterion re- 
quires the planet mass larger than 1.7 Mth with Hy = 
16H. This explains why a gap cannot be opened in our 
HD viscous runs (10% dip in VIO). However, a 50% deep 
gap is opened in our MHD run (M10B400) as discussed 
in §4.2, suggesting that the viscous criterion needs to be 
modified, at least for net vertical flux simulations. 

5.3. The Comparison with Previous results 

Our results are different from previous gap opening 
studies using net toroidal flux or zero flux MHD simu- 
lations (Winters et al. 2003, Nelson & Papaloizou 2003, 
Papaloizou et al. 2004, Uribe et al. 2011) in respect that 
previous simulations found similar gap depths between 
MHD cases and viscous cases 0. Furthermore, our main 
result that amax is higher in the whole gap region has not 
been reported before. Papaloizou et al. (2004) reported 
the increasing amax along the planetary wake in their 
Figure 19 & 21. But their Fig 21 suggests that a,nax 
decreases in the gap region. These differences between 
our results and previous results may suggest that differ- 
ent magnetic field geometries in the disk affect the gap 
opening process. 

To confirm that the deeper gaps in our simulations 
are due to the different field geometry we have applied, 
we have also carried out gap opening simulations with 
net toroidal magnetic fields. We have indeed observed 
similar planet-induced gap depths between net toroidal 
flux MHD cases and viscous cases, which is consistent 
with previous studies. 

More surprisingly, we found that, in net toroidal flux 
MHD simulations, amax decreases towards the gap re- 
gion (consistent with Figure 21 in Papaloizou et al. 
2004) , which is fundamentally different from the increas- 
ing amax towards the gap region in net vertical flux MHD 
simulations. This result will be presented in a later pa- 
per. 

* There are some subtle differences in these previous studies. 
Winters et al. (2003) found gaps are slightly shallower in MRI disks 
compared with viscous cases, while Nelson & Papaloizou (2003) 
suggest that gaps are deeper and wider in MRI disks. 



6. DISCUSSION 
6.1. Torque and Migration 

Planet migration is an important result of the planet- 
disk interaction theory (Goldreich & Tremaine 1980, Lin 
& Papaloizou 1993). As we have discussed in §4.1.3, in- 
teraction with a MHD turbulent disk changes the plane- 
tary torque density in comparison to HD cases. More im- 
portantly, it causes a "random walk" of the planet when 
fluctuations of the torque associated with turbulence 
dominates the Lindblad torque (Nelson & Papaloizou 
2004, Uribe et al. 2011). 

The time evolution of the total torques for runs 
M01B400, M03B400, and M10B400 are shown in Figure 
[HI The torques from both sides of the disk and the net 
torques are shown. The fluctuating torque is more im- 
portant for less massive planets (e.g. 0.1 Mt, M01B400) 
since the torque from the turbulence is proportional to 
the planet mass while the Lindblad torque is proportional 
to the square of the planet mass. 

For an intermediate mass planet (0.3 Mt, M03B400), 
the long lasting non-uniform background zonal flow (Jo- 
hansen et al. 2009, Simon et al. 2012) can also signif- 
icantly change the migration (Yang et al. 2009, 2012). 
This is shown in the middle panel of Figure [14] where 
the net torque is always negative due to the zonal flow, 
even though it should be zero if the background is uni- 
form in the shearing box limit. This net torque is almost 
20% of the total one-sided torque. However, caution has 
to be made by using shearing box simulations to study 
both stochastic migration and zonal fiows (Yang et al. 
2009,2012), and how zonal flow affects planet migration 
in global disks needs further study. 

It is unclear whether, in a very long time scale, the 
planet's stochastic motion can still lead to a one direction 
migration as in a laminar disk. We can indirectly address 
this question by averaging the one-sided torque over long 
times and comparing it with the Lindblad torque in HD 
disks. As shown in Figure [71 the averaged torques in 
net vertical flux MHD disks do not equal to the torque 
in HD disks. An excess torque very close to the planet 
is present around 0.2-0.3 H. A similar excess torque is 
found by Baruteau et al. (2011) (Figure 13) in global 
simulations. 

Several mechanisms can potentially explain this ex- 
cess torque, such as the "magnetic resonance" (Terquem 
2003, Fromang et al. 2005) and the horseshoe motion 
amplified corotation torque (Gullet et al. 2013). Unfor- 
tunately, we did not confirm any of these mechanisms, 
partly due to the complication of MHD turbulence. For 
the "magnetic resonance". We do not observe signifi- 
cantly density features associated with this resonance 
(Fig. 4 in Fromang et al. 2005). For the MHD coro- 
tation torque, our shearing box set-up do not exhibit the 
shift between the planet and the separatrices by design. 
This shift is global disks introduces the MHD corotation 
torque (Gullet et al. 2013). Thus, the mechanism behind 
this excess torque in our simulations is unclear. Exam- 
ining Fig. 3 in detail, we find that the circumplanetary 
region in MHD cases is asymmetric which may indicate 
some modification by the magnetic field to the circum- 
planetary region. 

6.2. Can we use a to study gap opening? 
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Although viscous disks have been extensively used to 
study gap opening and even planet accretion, our work 
suggests that caution has to be exercised when applying 
viscous models to real protoplanetary disks. 

First, a viscosity models MRI turbulence only in a 
statistical sense (Balbus & Papaloizou 1999). At any 
given instant in time, quantities such as the wake profile, 
gap shape, and torque can be very different in inviscid 
MHD turbulence compared to a viscous HD disk. This 
implies that if we are able to resolve the spiral structure 
induced by planets in protoplanetary disk with ALMA, 
the spiral arm may be distorted and disjointed which is 
very different from the smooth structure in viscous HD 
disks (Figlll). 

Second, even though the stress in MRI turbulent disks 
can be scaled to the disk pressure in a volume averaged 
sense, we find that since the turbulent stresses are not 
proportional to density, possibly due to global transport 
of magnetic fields, the value of the effective a can vary 
across the gap region. In our net vertical flux MHD sim- 
ulations, the gap region has a higher a than the rest 
of the disk. Calculating accurate models in viscous HD 
requires varying the value of a. An empirical formulae 
is presented in §4.2.1. Note that this formulae is de- 
rived from our simulations with limited net field strength 
range, planet mass range, and ignoring dissipation. Par- 
ticularly, the MHD turbulent strength can be affected 
by the dissipation coefficients (Lesur & Longaretti 2007, 
Longaretti & Lesur 2010) even in net vertical field sim- 
ulations. In order to understand the gap structure in 
the realistic simulation with dissipation, we need to first 
understand how magnetic fields are transported in such 
disks with density features. This demands numerical 
simulations with very high resolution, which is beyond 
the scope of this paper. 

Finally, properties of the circumplanetary region 
highly depend on the magnetic field geometry and field 
strength there, which is not isolated from the rest of 
the magnetic fields in MHD disks. Although the den- 
sity around the circumplanetary region is very high due 
to the gravity of the planet, the magnetic fields could be 
less affected by the presence of the planet due to the effi- 
cient global transport of the magnetic fields. If a ^ 1/2/3 
still stands there, the equivalent a in the cicumplanetary 
region could be very small. Moreover, the global geome- 
try of the magnetic field in this region means the stress 
is not completely local, but there are torques due to con- 
nections to the outer gap regions. It is unclear how this 
can be modeled in the local a formulation. 

6.3. Gaps in protoplanetary disks and observational 
Implications 

Current observations cannot constrain the structure of 
protoplanetary disks accurately due to the limited res- 
olution and disk's large optical depth at infrared. The- 
oretical calculations considering ionization network (e.g. 
Bai 2011) suggest the existence of the MRI "dead zone" 
(Gammie 1996). The stress (or a) can be very low in the 
"dead zone" making it similar to the inviscid HD disk. 
Based on our inviscid HD simulations, two gaps can be 
opened in the disk adjacent to a low mass planet. This 
"two gap" structure might be observable by ALMA in 
future. 

For MRI active regions in disks, the planet-induced gap 



structure depends on how the disk is threaded by global 
magnetic fields. Disks with net vertical fields have drawn 
greater attention recently since net vertical fields seem to 
be essential to produce disk winds and jets (Bai&Stone 
2012, Fromang et al. 2012, Lesur et al. 2012). If proto- 
planetary disks are threaded by net vertical fields, gaps 
can be opened by lower mass planets in contrast to the 
viscous HD models. 

On the other hand, if a gap is observed, "a models" 
may not be accurate enough to use the gap structure 
to determine the planet mass. Stratified MHD simula- 
tions with realistic disk ionization structures are needed. 
Eventually global MHD simulations are needed to de- 
rive a realistic gap shape. Some of these issues will be 
addressed in our next paper. 

Furthermore, the spiral arms produced by the planet 
can be fragmented by MRI turbulence. In a real pro- 
toplanetary disk, the spiral arms can be turbulent and 
dynamical, and different from the static spiral arms pro- 
duced in viscous "a models" . 

7. CONCLUSION 

We have studied the wake and gap opening by a low 
mass planet with both two dimensional HD and three 
dimensional MHD simulations in the unstratified local 
shearing box limit. Similar to the case of inviscid HD, 
the density wake launched by a planet steepens into a 
shock in the MRI turbulent disk. Thus shock dissipation 
dominates wave dissipation in MHD disks. Turbulent 
fluctuations cause the shock position to vary in time, and 
when averaged over hundreds of orbits the random wake 
positions form a smooth density profile which is remark- 
ably similar to that in viscous disks with the same total 
stress. On the other hand, we find the average torque 
density excited by the planet to be different from the 
torque density in HD simulations: it has a peak around 
0.2-0.3 H away from the planet compared to '~1 H in HD 
simulations. Furthermore, zonal flow makes the back- 
ground disk not uniform leading to a net torque even in 
the local shearing box limit. 

We have also studied gap opening criteria in both HD 
disks and MRI-driven MHD turbulent disks. In invis- 
cid HD disks, we find a 0.1 thermal mass planet can 
open two gaps at the position where the density waves 
shock, which subsequently deepen and merge into a sin- 
gle gap structure. This is consistent with the expectation 
from nonlinear semi-analytical calculations by Goodman 
& Rafikov (2001) and Rafikov (2002b), but contradicts 
the generally accepted "thermal criterion" for gap open- 
ing. For a disk with H/R=0.1, the timescale for this gap 
opening is around 10'^ orbits. In MHD turbulent disks, 
we find that a 0.3 Mth planet can open a 20% deep gap in 
a disk with equivalent a=0.08. Overall, both our HD and 
MHD results show that gaps can be opened by planets 
significantly below the thermal mass. 

By varying the strength of the initial vertical magnetic 
fields, we have studied the gap structure in MRI disks 
with different turbulent amplitudes, and compared the 
results to viscous HD disks with the same stress. The 
gaps in MRI disks having net vertical fields are deeper 
and wider in comparison to viscous HD disks. In partic- 
ular, we find that, with an initial vertical field /?o = 400 
in disks, the gap formed by a thermal mass planet con- 
tains of a density dip of 2, while only a 10% density 
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TABLE 1 

Hydrodynamic Models (2D) 



Case name Box size Planet mass Run Time Kinematic 



HxH Mt 27r/n Viscosity 

Inviscid DisJ<s 

n SxTB OT 400 

12 16x16 0.1 400 

13 32x16 0.1 400 
IlOinv 16x16 1 20 
lOSinv 16x16 0.3 20 
lOlinv 16x16 Ol 20 

Viscous Disks 

Vol 16x16 OT 100 0.113 

V03 16x16 0.3 100 0.113 

VIO 16x16 1 100 0.113 

VlOl 16x16 1 100 0.03 

VlOsv 16x16 1 100 0.056 

VlObsv 16x32 1 100 0.056 



dip is observed in the viscous HD disk. This differ- 
ence in gap shape is due to the fact that the MRI stress 
(r) across the gap depends only weakly on the density 
(AT/T — 0.25Ap/p in disks with initial vertical mag- 
netic fields) while the viscous stress is proportional to 
the density. 

Looking at this in a different way, the equivalent a 
parameter of the MHD turbulence is higher within the 



planet-induced gap of net vertical flux MHD disks. Re- 
markably, the correlation between the total stress mea- 
sured by a and the plasma /3 (~ l/2a) found generally in 
studies of the MRI (Hawley et al. 1995) still holds even 
inside the gap. Both turbulent magnetic fields and net 
vertical magnetic fields are more concentrated in the gap 
region with respect to the density, suggesting an efficient 
global transport of magnetic fields into the low density 
gap which leads to a variation across the gap. 

The fact that a is not constant within the gap, and 
that MHD stresses are not proportional to density across 
the gap (leading to deeper and wider gaps in MHD) calls 
into question the use of viscous HD simulations to study 
gap opening in protoplanetary disks. In the future, in or- 
der to model the dynamics of protoplanetary disks more 
accurately, we will study gap opening in stratified and 
layered MRI-turbulent disks. 
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Fig. 2. — The space-time plot for the disk surface density (y-direction averaged) with respect to x in the inviscid HD disk (case 12) having 
a 0.1 thermal mass planet at the center. X axis is in unit of orbits. Two gaps are gradually opened adjacent to the planet and they get 
deeper with time. The two dotted lines are where density waves excited by a 0.1 thermal mass planet become shocks. 



14 




Fig. 3.— The ^-direction averaged disk surface densities for MRI turbulent disks (left and middle panels, M10B400, M03B400, M01B400 
from top to bottom) and viscous disks (2D) having the same equivalent stress a=0.17 as the MRI disks (right panels, VIO, V03, VOl from 
top to bottom). The planets at the box center have different masses: 1 thermal mass in the upper panels, 0.3 thermal mass in the middle 
panels and 0.1 thermal mass in the lower panels. The left panels are the snapshots of the MHD simulations at a given time while the 
middle panels are from the same simulations but averaged over hundreds of orbits. Comparing the upper middle and right panels, we can 
see the gap is deeper in the MHD case. 
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Fig. 4. — The space-time plots of the j/2-direction averaged disk surface densities for M10B400, M03B400, M01B400, and B400 (top to 
bottom panels). A clear gap is opened in the M10B400, while density structures in lower mass cases (asymmetric with respect to x = 0) 
resemble zonal flow from B400 (the bottom panel). 
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Fig. 5. — The ^-direction averaged (Ax/c3)V ■ v of MRI disks with a 1 thermal mass planet at the center (upper left three panels, 
M10B400) and the MRI disk without planets (lower panels, B400). The leftmost panels show the color contours of (Aa;/cs)V ■ v at a given 
time. The second to left panels outline the contours with values smaller than -0.2, which identify the shock regions in disks. The second to 
right panels overlap these outlines in snapshots from 80 to 130 orbits with 5 orbits interval. The upper right panel shows the same outline 
for an inviscid HD disk (IlOinv). The planetary wake in the MHD case has similar shock strength as the wake in the inviscid HD case, thus 
in MRI disks the waves are dissipated by shocks. However, since the shock position varies with time, the averaged wake density profiles 
over hundreds of orbits are a lot smoother as shown in Figure \E\ 
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Fig. 6. — The (time and z-direction averaged) surface density profiles along y-direction at x=0.5,l,2,4 scale heights away from the planet 
in inviscid HD (lOSinv, dashed curves), viscous (V03, dotted curves) and MHD simulations (M03B400, solid curves). The soUd curves 

are scaled with the fast magnetosonic speed + v'^, detailed in §4.1.2). A remarkably good agreement is found between the averaged 

density profiles of MHD cases and those of viscous HD cases, even though the wave is dissipated by shocks in MHD cases (Fig. (5)1. 
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Fig. 7. — The time averaged torque densities excited by a 1 Mt, 0.3 Mt, 0.1 Mt planet in the disk (from left to right). The solid curves are 
from MHD simulations (M10B400, M03B400, M01B400). The dotted curves are from viscous HD runs (VIO, V03, VOl) while the dashed 
curves are from inviscid HD runs (IlOinv, I03inv, lOlinv). Torque densities in MRI disks have peaks closer to the planet. 
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Fig. 8. — The j/z-dircction and time averaged disk surface densities (upper left panel), Maxwell Stress (upper right panel), equivalent 
Q-Max (lower left panel), and net vertical flux (lower right panel) for vertical net flux MHD simulations with a 1 thermal mass planet at 
the box center. Regions closing to the planet {\y\ < H) arc masked out for the averaging. \x\ < H is shown in the upper right and lower 
panels as the shaded region. Two different initial field strengths have been applied (M10B400(dark curves) and J\110B1600 (light curves)). 
The equivalent viscous cases(V10sv, VlOst) are over-plotted as the dotted curves. In order to be compared with the Maxwell stress and 
OMax, the viscous stress (Txy) and a are multiplied by 3/4 in these plots. Clearly, the gaps are significantly deeper in MHD cases than 
the viscous cases since in MHD cases the stress is flat compared to the density, or, in other words, a increases towards the gap region. In 
the lower right panel, the initial vertical field strengths axe plotted as the flat solid lines. Net magnetic fields are concentrated in the gap 
region. 
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Fig. 9. — The disk surface density (left panels), and the ratio between the stress and density {a, right panels) for viscous (VIO, upper 
panels) and MHD (M10B400, lower panels) cases. For the MHD case, xy component of the stress is the Maxwell stress. For the viscous case, 
the viscous stress {Txy) is multiplied by 3/4 to be compared with the Maxwell stress. Compared with the viscous case, the stress-density 
ratio (a) in the MHD case is higher in the gap. 
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Fig. 10. — Similar to Figure [8] but for more MHD and viscous cases. Upper panels: the case with a more massive planet (M30B400), 
and a y-direction longer boxes (M10B1600b)). The viscous case (VlObsv) is shown as the dotted curve. Bottom panels: the cases with a 
less massive planet and various net vertical field strengths {M03B400, M03B1600). The dotted curve is the viscous case V03. In all these 
cases, a of MRI disks peaks towards the gap region. 
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Fig. 11. — Left Panel: Volumetric rendering of the density during the interaction between the planet and MRI turbulent disks for 
MlOBlGOOb. Regions with densities below 0.8 midplane density are removed. Right Panel: The time averaged aMax at the z=0 slice for 
this case. The time averaged magnetic field streamlines are also plotted. The turbulent component of the magnetic fields is averaged out 
and the averaged fields show the net vertical field geometry. The net vertical magnetic fields in M10B1600 seem to diffuse freely into the 
gap, causing a higher a. 
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Fig. 12. — Left panel: the excess of the ?/2:-direction and time averaged Maxwell stress with respect to the excess of the averaged density 
at each x position across the gap region for the M10B400 (black curves), M10B1600 (red curves), M30B400 (blue curves), and M10B1600b 
(green curves). < B^By >o is the Maxwell stress at the position where the averaged p = po = I. The dotted curve from the lower left to 
the upper right shows the stress-density relation in a viscous disk AT/T = Ap/p, while the best fit for the stress across the gap in MRI 
disks is AT/T oc 0.25Ap/p. Right panel: The j/z-direction and time averaged /3 with respect to the averaged disk Maxwell stress aji/aa: at 
each X position across the simulation domain. The point moves to the lower right when it crosses the gap where it has lower /3 and higher 
'^Max- The solid line is /3 = l/2a. 
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Fig. 13. — The space time plot for the j/z-direction averaged disk surface densities for the M10B1600 (net vertical flux), and VlOla 
(viscous HD). However, the planet's potential is suddenly changed to be zero in the middle of the simulation. As shown, the gap feature in 
M10B1600 persists longer than that in a viscous disk. This is consistent with the fact that the stress across the gap in MRI disks is more 
uniform than the density, so that any existing density feature takes a longer time to be altered by turbulent stress in net vertical flux MHD 
disks. 




Fig. 14.— The integrated torques with time for M01B400 (1 Mth), M03B400 (0.3 Mth), M10B400 (0.1 Mth) (left to right panels). The 
dashed curves are the torques from x > side of the disk while the dotted curves are from x < side of the disk. The solid curves are 
the net torque. These curves are averaged over 1 orbit period. The red curves are the net torque averaged over 10 orbits period. For 
less massive planets, the torque due to the turbulence can dominate the Lindblad torque, leading to the planet's random walk. For the 
intermediate mass planet (0.3 Mt case), the net torque is actually non-zero (negative in this case), contradictory to the expectation that, 
due to the symmetry of the shearing box, the averaged net torque should be zero. This non-zero net torque is caused by the non-uniform 
disk background from the large scale MRI zonal flow. 



